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Due to recent successes of a statistical-based nonlocal continuum crystal plasticity theory for 
single-glide in explaining various aspects such as dislocation patterning and size-dependent plasticity, 
several attempts have been made to extend the theory to describe crystals with multiple slip systems 
using ad-hoc assumptions. We present here a mesoscale continuum theory of plasticity for multiple 
slip systems of parallel edge dislocations. We begin by constructing the Bogolyubov Born-Green- 
Yvon-Kirkwood (BBGYK) integral equations relating different orders of dislocation correlation 
functions in a grand canonical ensemble. Approximate pair correlation functions are obtained for 
single-slip systems with two types of dislocations and, subsequently, for general multiple-slip systems 
of both charges. The effect of the correlations manifests itself in the form of an entropic force in 
addition to the external stress and the self-consistent internal stress. Comparisons with a previous 
multiple-slip theory based on phenomenological considerations shall be discussed. 
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I. INTRODUCTION 

Statistical mechanics provides an optimal framework 
and various tools for studying emergent phenomena 
from a complex conglomerate of bodies — may they be 
molecules of gases, polymer chains of rubber, or crys- 
talline defects. The use of correlation functions in 
analysing two-dimensional solids and their defects has 
been proven very successful in the past. For example, 
Mcrmin showed that two-dimensional crystals do not 
have conventional long-range order, but can have "di- 
rectional long-range order." ^ Nelson et al. applied the 
technique to explain dislocation-assisted melting in two 
dimensions. ^'^ Over a decade ago, Groma proposed a 
theory to describe dislocations and their motions using 
distribution functions and probability arguments."^ Un- 
like the existing continuum theories at the time,^*^ the 
new formalism was physically motivated and incorpo- 
rated correctly the long-range nature of dislocation in- 
teractions. Several variations of this work — all of which 
reduce to the same two-dimensional theory — also exist 
for three dimensional dislocation systems. 

Although having laid out the foundation for possible 
interactions of many-dislocation configurations, Groma's 
pioneering work did not investigate these correlated ef- 
fects in details. Zaiser et al. considered explicitly 
the evolution of dislocation correlations by extending 
Groma's theory for systems of single-slip, parallel edge 
dislocations.^ They were able to qualitatively obtain the 
correct scaling behavior of the evolution equations for 
both single and pair correlation densities, and explained 
some general properties of these functions. Their formu- 
lation, however, was limited to only one active slip sys- 
tem and the analytical forms of pair correlation functions 
were not derived. In a later work, Groma et al. included 
the influence of dislocation correlations in the form of a 



local back stress. '^'^ Yefimov et al. connected this statis- 
tical description to a continuum crystal plasticity theory 
and applied this to various boundary value problems. ^^'^^ 
While the theory successfully captured most features ob- 
served in discrete dislocation simulations, its ad-hoc ex- 
tension to multiple slip systems failed to explain size ef- 
fects in single crystal thin films. The main goals of this 
paper are: (1) to correctly describe and obtain analyti- 
cal expressions for dislocation pair correlations, and (2) 
to systematically generalize the approach of Groma et al. 
to multiple slip systems. 

We begin, in Sec. II, by introducing ensembles of dis- 
locations and deriving the partition function for mul- 
tiple slip systems. The n*'^-order dislocation densities 
and dislocation correlation functions are subsequently de- 
fined. We construct the Bogolyubov-Born-Green-Yvon- 
Kirkwood (BBGYK) integral equations in Sec. III. These 
equations link correlation functions of order n to those 
of order n + 1 (a technique generally used in the study 
of dense gases and fluids). The integral equations are 
expanded in powers of interaction strength (the ratio be- 
tween the interaction energy and 'thermal' energy). We 
then obtain a set of approximate integral equations for 
pair (n = 2) correlation functions after applying a closure 
approximation to truncate the series. These equations 
are valid regardless of the form of the interaction poten- 
tial, and thus are applicable to other systems, provided 
that this pair interaction vanishes at a large distance. 

By appealing to Peach-Koehler interaction, analytical 
expressions for pair dislocation densities for single and 
multiple shp systems are derived in Sec. IV and Sec. V 
respectively. Our single-slip solution agrees with the re- 
sult from the study of induced geometrically necessary 
dislocations (GND) in terms of a single pinned disloca- 
tion by means of a variational approach.^'* The disloca- 
tion spacing l/^J~p emerges as a natural lengthscale in 
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this formulation in accordance with the scahng study by 
Zaiser et al.'' Our analysis further shows long-range at- 
tractive correlations when more than one slip system are 
present, confirming the absence of dislocation patterning 
in single glide systems as observed in many discrete dis- 
location simulations^^"^^ and explained in a recent three- 
dimensional continuum plasticity theory.^'^ 

In Sec. VI, we write down the transport equations for 
both total disloc;ation densities and GND densities on 
each slip system under the influence of Peach-Koehler 
forces from both single and pair dislocation correlations. 
While the former gives a self-consistent, long-range in- 
ternal stress contribution, the latter exerts an additional 
short-range, entropic force due to a deviation away from 
a preferred dislocation arrangement in the form of a back 
stress. The formulation is a direct extension of the work 
by Groma and Zaiser**'^'^" for crystals with one active slip 
system. Using knowledge of the pair correlation func- 
tions, we obtain a complete description of the back stress 
as a function of slip orientations — which previously had 
been incorporated using ad-hoc phenomenological con- 
siderations in the multiple-slip theory. ^^'^^ 

Finally in Sec. VII, we contrast our theory with the 
multiple-slip theory of Yefimov et al.^^'^'^ While both the- 
ories propose that interactions among slip systems de- 
pend solely on relative angles of slip orientations, the 
functional forms are different. We attribute the failure 
of the earlier theory in explaining size effects in single 
crystal thin films partly to this difference and partly to 
the treatment of dislocation nucleation in the theory. 



II. DEFINITIONS OF THE BASIC QUANTITIES 

Consider a system containing r species of disloca- 
tions and denote the coordinate of the dislocation 
of species s by is- The dislocation configuration {N} 
is the set of the coordinates of all dislocations, where 
N = {Ni,N2; N3, N4;...; Nr-i,Nr) denotes the "collec- 
tion" of dislocations of type s. In this convention, odd 
and even slots respectively contain plus and minus dis- 
locations on distinct slip systems. We introduce the 
notation {N -|- Ig} to denote the addition of an extra 
dislocation of species s to {N}, while similarly a config- 
uration {N} with coordinates of n removed is indicated 
by {N - n}. 

The interacting Hamiltonian of the system can be 
written as the sum of potentials ^(isi — iss) of all pairs 
of dislocations 

c^n({n})= j2 (1) 

Sl<S2 i<j 

We can define a canonical partition function of configu- 
ration N by 



-!7N/feBT 



d{N}, 



(2) 



where the integrations are taken over the "volume" mea- 
sure d{N} = Yrg^^(flsd'^2s---(PNs of the dislocation 
configuration at {N}. Consider the coordinates of a par- 
ticular set {n}, the probability of observing the configu- 
ration n in d{n} about the points in {n} irrespective of 
the remaining coUecttion N — n is 



P(f)({n})d{n} 



d{n} 



-C/N({N})/fcBT 



d{N - n} , 



(3) 

where / P^\{n}) d{n} = 1. The probability density 

of observing any statistically equivalent possible collec- 
tion n within the volumes d{n} about the points {n} is 
therefore 



P^."^({n}) = n 



ATI 



LI (AT, - n«)! 



^',^"^({n}). 



(4) 



By using Boltzmann distribution, we assume that our 
system is ergodic. and thermal equilibrium exists and can 
be reached. System of dislocations which drift along the 
local force (thus implying glide to be accompanied by 
some amount of climb) subject to thermal noise would 
certainly fit the criterion. 

Consider now an open system (which could be real- 
ized, say, by allowing for nucleation and annihilation of 
dislocations as the system relaxes); a grand canonical 
partition function is given by 



-=En 



N>Os=l 



NJ 



--N : 



(5) 



where Zg is the activity of species s. The prefactor arises 
from integrating away the nionientuni degrees of freedom 
in the Hamiltonian which are irrelevant to this problem. 
The probability V of the occurence of configuration N in 
the open system is therefore 



s=l 



NJ E 



(6) 



Finally, the probability density of observing any ni dis- 
locations of species 1, n2 dislocations of species 2, etc., 
{any collection n) in d{n} at {n} is 



pW({n})= ^T'npJ^HW). 



(7) 



N>n 



The summation is taken over all collections N greater 
than or equal to n, i.e., for all N\ > m, N2 > n2, etc. 
We take Eq. (7) as the definition of dislocation density 
of order (n). Explicitly we have 



N>n s=l ^ * 

'■e-^N(W)/teT^{N-n} (8) 
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This definition of an (n)**^-order dislocation density is 
equivalent to the ones used by Groma"' and Zaiser^ in 
the realization of an open system. Finally, we define 
the (n)*^-order correlation function g^'^\{n}) through 



(W) 



s=l 



5^"H{n})- 
(9) 



III. 



DERIVATION OF THE BBGYK INTEGRAL 
EQUATIONS 



The Bogolyubov-Born-Green-Yvon-Kirkwood inte- 
gral equations first appeared in the study of classical flu- 
ids with a total potential energy given by the sum of pair 
interactions. ^^"^^ They provide a set of relations between 
distribution functions of fluid density at different orders. 
Here we extend the BBGYK formalism to include the 
non-central interactions of dislocations in a multicompo- 
nent system. ^^■^'^ We proceed in three steps: (1) take a 
derivative of the (n)*'^-order dislocation density with re- 
spect to the coordinate of one particle of the interested 
species; (2) express the result in terms of the next higher 
order densities; and (3) convert the integral equations of 
densities into those of correlation functions. 

Differentiating p(")({n}) as expressed in Eq. (8) with 
respect to dislocation 1 of species 1 located at li we find 



where we absorb 1 /k^T into the definition J7n := U/k^T. 
The derivative of the potential can be separated into two 
parts: 

r ris r 

8=1 1=1 8=1 i=ns-\-l 

{i,s)7t{l,l) 

(11) 

Direct substitution of Eq. (11) into the integrand of 
Eq. (10) splits the expression into two integrals Ii and 
l2- Notice in the first integral that the derivative of the 
potential does not depend on coordinates {N — n}, and 
thus can be taken out of the integral, yielding 



h = -p("n{n}) E E - (12) 

s=l i=l 
(i,s)7^(l,l) 



with the aid of Eq. (8). The second integral I2 requires 
a little more work: 



r r 



AT,, 



s=l N>n 



f Vr^w(li-i'.)e-^-«^»4N-n} (13) 

^ A— nr. I 1 



i=ns+l 



vry"H{n})=-^E[ri 



N>n s=l ^ * ^' 

I e-0^(i^})v^U^{{-N}) d{N - n} , (10) 



The expression involves integrating ig over the sample 

size. Since each integral over ig between + I < i < Ng 
is equivalent in infinite space, the summation therefore 
gives a factor of (Ng — Ug). The remaining integrals over 
all other dislocation coordinates are unaffected. 



Eq. (13) thus becomes 



r r 

^^=-EiE[n 



2=' 



s=l N>n s' = l ^ * ^ ' 

J V^uiU - K + lj.){ J e-^-(W) d[{N - n}\{(n, + ljj] } d^K + 1). 



= -E/vr.^(i-K + iw{^E[n^ 

8=1"' '-"N>n s'=l ^ 



(iVs ~ Ug) 



ng,)\ 



E/vr.^ 

s=l •' 



J d[{N - n}\{(n, + l) j] | d\n, + l) 



u{U - {ng + l)g) p("+i=^({n + 1 J) d'ing + l)g 



(14) 



The symbol d {N — n}\{(ns-|-l)s} represents the volume measure of {N — n} without d {ns + l)g. Collecting both 
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Ii and I2 from Eq. (12) and Eq. (14), we arrive at the BBGYK equations for the (n) -order dislocation density: 

Vr^p(")({n}) = -p("H{n}) ^£^Vr^u(ii -4) / Vr^u(ii -(;^.)p("+i'H{n + iJ)d'K+i)« (15) 

One can obtain a series of integro-difFerential equations for the correlation functions g*-"^ from Eq. (15) by expanding 
out p^"^({n}) using Eq. (9). All but two of the single dislocation densities on the left and right-hand sides of the 
equality cancel which results in 



{r,ns ) 

p(li)5("n{n})] = -p(li)ff("H{n}) E ^1,^(11-?.) 



r „ 

p(li)^ J %uil, - in, + l),) piin, + l),)g^'^+^^\{n+ 1,}) d^{n, + l), . (16) 



The first order densities p(l) that plague the expression 
can be removed by first using the product rule to the 
left-hand side (LHS), then dividing both sides by p(l). 
The LHS becomes 



LHS = Vr^g(")({n}) + ff(")({n}) 



, Vr,p(ri) 



The ratio of the derivative of the first-order density with 
itself can be rewritten using Eq. (15) specialized to first 
order, giving 

where ^ = {ns + l)s is the position of the (ng + iy^ 
dislocation of species s, and g^'^\li,^s) represents the 
pair correlation function between the first dislocation of 
species 1 at li and the (n^ + 1)"^ dislocation of species s 
at ^g. This expression could be incorporated seamlessly 
into the right-hand side of Eq. (16). The final result is^^ 

(r,ns) 
(s,i)7^(l,l) 

-E / vr,H(ii-i;)p(i;) 

gi-+^^){{n + l,})-g(^\{n})g^'\iui)] d% . 

(17) 

For the remainder of this paper, we shall restrict our 
attention to the Peach-Koehler interaction. Recall that 

the interaction energy between two parallel edge disloca- 
tions of length L (over thermal energy k^T) in an infinite 
medium is^'^ 



where T 



27r(l - v)kBT 



and 



{fit ■ {is - Js')j (mi, ■ (is - Js')) - 



+ 



■ (19) 



Here denotes the slip-plane normal of species s. The 

relative strength T represents the ratio between disloca- 
tion interaction energy versus the system's thermal en- 
ergy. Note that the latter originates from the use of 
Boltzmann distribution in Eq. (2) to describe the equilib- 
rium configuration of systems with thermal noise. It was 
pointed out by Groma et al.^"* that, in systems where 
dislocations are confined to their slip planes, the glide 
constraint acts as an effective temperature preventing 
the systems to relax by means of (lisloc;ation annihila- 
tion. Seen in this light, the temperature T in this theory 
is not physical temperature but a fictive temperature as- 
sociated with the disorder in dislocation distributions.''''* 
As the dislocation configuration becomes more and more 
correlated, F becomes smaller. 

For an explicit dependence on F to use as an expansion 
coefficient, we rescale the distance by the square-root of 
the relative strength, vT'r 1— > r. Eq. (17) specialized to 
second order gives 



Vj-ff(2)(i,2)=F5(2)(i,2)Vj-^(l,2) 

+ E / Vr^(l,3,)p(3,) 



g^^\l,2X)-9^'\W\lX) 



d% . (20) 



u{is - 3s') = -F 'il}{is - js>) 



(18) 



Here we have simplified the notation even further by sup- 
pressing all irrelevant subscripts: vectors 1 and 2 simply 
denote the positions of dislocations 1 and 2 with their 
corresponding species. The summation is taken over 
all s species present in the system. 
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Wc proceed by assuming that the correlation functions 
have the foUowing forms: 



5(2)(i,2) = l + r/(2)(l,2), 



(21a) 



5(3)(1,2,3) = i + r 



/(2)(l,2) + /(2)(l,3)+/(2)(2,3)' 



+ rV('^(i,2,3), 



(21b) 



for any vectors 1, 2, and 3. The functions /'^^(l, 2) 
and /'^'''(l, 2, 3) should asympotically vanish along the 
boundaries of the sample, or as |1 — 2|, |1 — 3|, |2 — 3| — > oo 
for an infinite system. Note in particular that 

^(3) (1, 2, 3) - (1, 2)5(2) (1, 3) = r (2, 3) 

+ r2[/(3)(l,2,3) -/(2)(r,2)/(2)(r,3)l. (22) 



So far no approximation has been made. The Eqs. (21) 
governing the second-order correlations naturally involve 

the third order correlations. To systematically close the 
chain at the second order, we substitute Eqs. (21) and 
(22) into Eq. (20) to produce a set of integro-differential 
equations of /(^^ and /(■^'' for each power of T. This 
technique was introduced by Bogolyubov in the study of 
correlations in Coulomb interactions^^ and has since been 
widely used in both high energy and condensed matter 
communities in renormalization group theory. 

The equation of power r° gives an identity. After in- 
tegrating away Vj because /(^^ and ijj vanish along a 
boundary, the equation of power T becomes, 



/(2)(1,2) = ^(1,2) + ^ / V'(l,3,)p(3,)/(2)(2,3,)d% 

s=l 

(23) 

This equation is the key result of the analysis. In the 
following sections, we shall use it to obtain dislocation 
pair correlation functions for systems with one (Sec. IV) 
and many (Sec. V) active slip systems. 



IV. PAIR CORRELATION FUNCTIONS FOR 
SINGLE SLIP 



To illustrate the use of Eq. (23), we first apply it to the 
case of one slip system containing two species of disloca- 
tions (denoted + and — ). According to Eq. 19 valid for 
an infinite sample, ■i/'(l, 2) = "ipil — 2) — ip{2 — 1) which 
implies that /"^(1,2) = /"^(l-2) = /"''(2-1). Without 
loss of generality, we can take the origin to be at 2 and 
thus, from (23), we obtain the following set of integral 



equations: 

r\f) = ^i(f=) + j (ff'Mr- f') 

[p\r')r{r')-p{r')r-{r')] 
/-(r) = -Vi(r1 - y d2f>i(f - f') 



[p-(f ')/--(r") - P*{r')r{f') 



r (r) = + j d^r' V'l {f-f') 

[p-{f')r{r') - p\nr{r')] 

f-\f) = -^i(r) - J (ff'Mr-r') 

[p\r')r{r') - p-innn] 

In the current context, Eq. (19) reduces to 



^i(r) = V^^(r) = -V^-(r) = ln(|r=1) + 



\f\ 



2 ' 



(24a) 
(24b) 
(24c) 
(24d) 

(25) 



where we orient our {x,y) coordinate system in such a 
way that the slip direction points along the x direction. 
The minus signs in Eq. (24) arise from a sign difference 
in the interactions between plus-plus dislocations versus 
plus-minus dislocations as shown in Eq. (25). By com- 
paring Eq. (24a) against (24d), and Eq. (24b) against 
(24c), we find that f*^{^ = -f'^ir) and f*-{r) = 
—f~~{r). These symmetries further imply that f*'^{r) 
■ Finally we obtain 

f*{f) = V'i(r) 

+ j Mr- nnn [p\f') + p-(f o] <er' . (26) 

Our general formulation in the previous section allows 

for spatial variation of an uncorrelated density p{rg). 
Without externally applied force, p{rs) — (Ng) /A is con- 
stant in space. An analytical solution to Eq. 26 can be 
obtained for constant p* and p'. The dimensionless na- 
ture of the interaction potential ipi suggests a change of 
variable ^/p"^ + p~ r ^ 'r (note that p* and p~ are always 
positive). The resulting dimensionless integral equation 

r^(r) = ^i(r) + j Mr- nr{r')d^r' (27) 

can be solved directly by applying = (9^ -|- dyY on 
both sides of the equation and using the identity 

AVi(r) = 2-K/:^5{r)+2'K{dl-dl)6{r) = 47ra^5(r). (28) 

Eq. (27) then becomes 

= 47r92[/- + <5(r)] , (29) 
whose explicit solution is 

y 

f** = - smh.{y/T^y)Ki{y/TTr) — cosh(V7rj/)/sro(-\/7rr) , 
r 

(30) 
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with Kq{-) and Ki(-) the zcroth and first order modi- 
fied Bessel functions of the second kind. With the aid 



of Eq. (21a), the correlation functions g 



(++) 



original coordinates, 



9' 



) and 



*\ correct to 0{r^), can be expressed in the 



(r) = 1 + r 



(r) = l-r 



y 



smh{koy)Ki{kQr) 



cosh.{koy)Ko{kor) 



siioh{koy)Ki{kor) 



cosii{koy)Ko{kor) 



(31b) 



where ko = \/t^T{p^ + p') gives an inverse "Debye ra- 
dius" of the dislocation cloud. The third order corre- 
lation functions correct up to O(r^) follow straightfor- 
wardly from Eq. (21b). The validity of Eq. (31) can be 
verified by comparing g'^**'> (r) — g^*'^ (r) with the disloca- 
tion difference, or GND, field K(r) in Eq. (15) of Ref. 14. 
In this latter work the same expression is obtained for the 
induced GND due to a single pinned dislocation, which 
was interpreted by the authors as the pair correlation of 
dislocations in a relaxed system. 

It is interesting to note that the pair correlation func- 
tions depend only on the scaled space coordinate y/pf 
{p = p* + p~ being the total dislocation density) in agree- 
ment with the scaling argument given by Zaiser et al.® 
This dependence also holds in the multiple-slip case to 
be discussed in the next section. 



By direct substitution of + and — into a and &, it is 
immediate that f]*{r) = -fi]{r) and /T:(r) = -fTj{r), 
which further implies that 



N 



k=l 



i'ik * [Pkfjk + Pkfjk] ■ (33) 



(31a) With this, Eq. (32) reduces to 



fij = Aj + '^i^ik* [Pkfjk] : 



(34) 



fc=l 



where the superscripts have been omitted and p^. = 
Pk + P~k total dislocation density of both types on 

slip k. We thus effectively reduce the number of coupled 
equations to N"^ . Note also that because of = ipji, 
there are only N{N + l)/2 independent t/jij's. 

As seen from the single-slip c;ase, Eq. (34) subjected to 
an arbitrary distribution of the local density pk{r) can- 
not be solved analytically. For spatially independent pk, 
however, these equations can be decoupled. Let be the 
relative population of density in slip system k relative to 
the total density p, i.e., pk = XkP where Ylk=i = 1- 
We can then perform a change of variable r r to 
absorb the p-dependence. In addition, in Fourier space 
(indicated by a superposed '~), a convolution becomes a 
product. We can solve the Fourier-transform of (34) for 
fij by essentially performing a matrix inversion on 



(35) 



V. PAIR CORRELATION FUNCTIONS FOR 
MULTIPLE SLIP 

The procedure to obtain the correlation functions for 
a system with multiple slips follows the same types of ar- 
guments and expansions as those for single slip. We shall 
further develop the integral equation (23) for a system of 
A'^ slip systems, each with two charges, and subsequently 
give an explicit analytical solution for the pair correlation 
function in the case where the difference in slip orienta- 
tion angle between adjacent slip planes is constant. 

For an A'-slip system with both types of charges, we 
have 4A''^ coupled integral equations for different pairs 
of 1 and 2 in Eq. (23). To reduce the number of equa- 
tions, and essentially decouple them, some symmetry ar- 
guments can be employed. For an infinite system, 

^ = ^^-^ = -^1^^ and i>t^r^, 

where the superscripts denote the charges of the first and 
second dislocations, while the subscripts show the slip 
systems in which they live. Eq. (23) can be re-cast using 
the convolution operator * and the symmetry of V't^ as 



N 



(32) 



The Fourier representation of ipij in Eq. (19) can be 
expressed very simply in polar coordinates {k,(pk), 



47r 



sin(0fc - 9i) sin(^fc - 9j) 



"F 



{fhi ■ k) {rhj ■ k) 
(36) 

where 6i is the angle that slip plane i makes with the x 

axis (which can be chosen arbitrarily, so that Oi = iir/N). 
Owing to the simple form of (36), the solution to (35) is^^ 



(37) 



k=\ 



where we have used X]„ i^in^Pn] = i'lj I]„ Am- Eq. (37) 
shall be used in the derivation of the evolution law for 
parallel edge dislocations in a multislip system in the next 
section. 

To verify that Eq. (37) is appUcable in glide-controlled 
systems, we consider an ensemble of 1500 relaxed config- 
urations of 64 plus and 64 minus dislocations randomly 
placed on a 1 /zm^ square and restricted to move along 
their glide directions. The simulations were performed 
with periodic boundary conditions in the absence of ther- 
mal noise. The glide constraint helps prevent dislocation 
annihilation, and thus, to fix the total number of disloca- 
tions and to maintain the finite effective temperature. As 
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(a) 



-3 -2 -1 

(b) 



o|ll 



i 



FIG. 1: (a) Discrete dislocation result and (b) theoretical 
prediction of the correlation function /12 between plus dis- 
locations on 60° and 120° slip systems. Values increase to- 
wards brighter regions. Coordinates are measured in units of 
1 / yp. Dashed lines indicate the two slip directions where the 
plus-plus anti-correlation is underpredicted due to the glide 
constraint of the discrete dislocation simulations. The fit- 
ting parameter due to rescalings of length was found to be 
ko ~ 22^. 



an example, Fig. 1 shows (b) the density plot of the theo- 
retical correlation function between plus dislocations 
on 60° and 120° slip systems against (a) the simulation 
result. The erroneous oscillations in Fig. 1(b) along 0° 
and 90° lines are caused by the numerical inverse Fourier 
transform operation of Eq. (37) . (The general closed form 
solution of a double-slip pair correlation function does 
not exist for an arbitrary pair of slip orientation angles.) 
Overall, the theory gives accurate angular predictions ex- 



cept along the two slip directions where it underpredicts 
the same-sign anti-correlation due to the suppression of 
climb. The plot of the correlation function along the st 
axis is shown in Fig. 2. Very close to the origin, the 
function diverges logarithmically as does the unscreened 
potential. About one dislocation spacing from the core, 
the correlation function decays as 1/x'^. 




1.0 1.5 2.0 2.5 3.0 ^ ' ^ 



FIG. 2: Cross-sectional plot of the data points versus theo- 
retical curves of the pair correlation function (Fig. 1) along x 
axis. After a short distance away from the core, the function 
has a power law decay of as shown with the dashed line 
in the log-log plot in the inset. 



The real-space solution to Eq. (37) is possible if we 
assume that the angle between each adjacent pair of slip 
planes is constant. For any N G Z+ and > 1, 



N 



ra=l 



~N~) 



N 
2" 



(38) 



regardless of (j)k. With the above identity, the denom- 
inator of fij becomes angular independent and can be 
integrated directly. The final result, with 



reads 



ko=^2TTNTp, 
-2 r cos(20 -9,- 9j) cos{9, - 9j) 



kor 



- sin(0 - 9,) sin(0 - 9j) K2ikor)j. (39) 

At large distances, the first term dominates and the 
pair correlation decays like 1 /r^ (except along the direc- 
tions where the argument of the cosine is 7r/2, 37r/2, etc.). 
Compared to the single slip case (Eq. 30) where the pair 
correlation diminishes exponentially (except along the 
dislocation wall direction), the presence of extra slip(s) 
suppresses the Debye screening. It should be noted that 
—fij{r) can be thought of as the effective interaction po- 
tential due to screening. More precisely, ^ ^ fij 
is the Peach-Koehler force felt by a positive dislocation 
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on slip system i due to the induced screening of dislo- 
cations on slip system j. It has been shown^^ that, for 
single-slip system, the attractive parabolic potential in 
the glide direction (taken to be along x) falls off with a 
prefactor of along the wall direction. Series ex- 

pansion of (p in Eq. (39) about di and 9j reveals that, for 
multiple-slip system, the prefactor of the paraboHc po- 
tential about the glide directions decays as 1 /r^ — slightly 
more slowly than the single-slip case. This could explain 
the necessity to include more than one slip system to see 
the formation of cell walls and grain boundaries in two- 
dimensional discrete dislocation simulations prohibiting 
climb motion. ^'"'"^^ The analysis also confirms the "direc- 
tional long-range order" of two-dimensional crystals as 
rigorously proven by Mermin.^ 



VI. DERIVATION OF A MULTIPLE-SLIP 
EVOLUTION LAW 

To arrive at a set of transport equations for an ensem- 
ble of multiple-slip dislocation systems, we extend the 
treatments of Groma et al. in Ref. 4, 32, and 9. The evo- 
lution eriuations for the uncorrelated single-dislocation 
densities on slip system i read: 



dtp;{fi,t) = -{bi-^) 



ind 
ij 



(40a) 



+ 



(40b) 



where the dislocation mobility has been absorbed into 
the rescaling of time t. With the assumption that all dis- 
locations have the same magnitude b, the Burgers vector 
can be written as bi = 6s, (sj and m, respectively are the 
slip direction and slip plane normal direction of slip sys- 
tem i). T^"'^(ri — Vj) is the resolved shear stress exerted 
on a dislocation at on slip z by a dislocation at on 
slip j, and can be written as 

Tlf{r} = Si-crymi = Gb {si ■ V) (m^ • V) {rrij ■ V) [r^ In r] . 

(41) 

Here, G = /i/(27r(l - v)) = £;/(47r(l - v"^)), where E, /i, 
V are the Young's modulus, shear modulus, and Poisson 
ratio respectively. 

Addition and subtraction of Eqs. (40a) and (40b) give 
the evolution equations for the total dislocation density 



pi = p'l + p • and the GND density Ki = p* - p^■. 



dtPi = -{h-V) 



j 



+ V. / d^rj {pI] + p1] - pI] - Pij ) 



ind 



(42a) 



9tKi = -{bi- V) 



PiTi 



+ 



Yl / '^'^J (PiJ - - Pi] + Pi]) 



ind 



(42b) 



In accordance with (9), the dislocation-dislocation den- 
sity can be written as 



pTj = PKri}pj {rj)gtj in - fj) 

= Pi(rl)p/(rS)(l + C'(rl-fS-)), 



(43) 



where s,s' G {+)— } and, according to (21a), dfj = 

r/f| . In terms of the single and pair correlation func- 
tions, the total dislocation density and GND are 



Pif =Pi{n)pj{rj) + \{- Pi{ri)pj{fj)d^j 
+ Pi{ri)Kj{fj)[cl^j + dtj] 

+ Ki{fi)pj{fj)[d^j - d^j] + Ki{fi)Kj{fj)d^jY 



(44a) 



4f = i^i{n)i^j{rj) + ]^[pi{ri)pj{rj)[(Pi^ - df^-] 

+ Pi{n)Kj{rj)d1^ - Ki{ri)pj{rj)d1^ (44b) 
+ i^i{ri)Kj{riM^+dl^]], 



where (if - 



(l/2)(4: + d:;), and = 
(1/2)(4t - dTj-). After substitution of Eqs. (43)-(44), 
Eq. (42) becomes 



dtPi = -{bi ■ V) [Kiirr' + - r| - r^) + p.rf] 



(45a) 



dtKi = -{bi ■ V) [piirr + rf - rl - r^) + K^rf] 



(45b) 
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in which 



(46a) 



(46b) 



(46c) 



< = I E / - - ri) 

^ J J (46d) 

The term ^ = (1/4) (rfj-} + rf"; + + rf-;) in (46b) 
involves averaging over pairs of correlation functions. 

Terms involving rf in Eq. (45) can be cast away by 
going into a "co-moving" frame of pi and Ki respectively. 
Although = d~~ = = ^d~ij and hence d\j should 

vanish by definition, this is hardly the case when, e.g., 
the system is strained through external loading. Only 
one of these correlation functions dominates locally, re- 
sulting in a nonzero d\y Similarly the contribution from 
flow stress, r/, is greatest in regions with equal popu- 
lation of plus and minus dislocations; in most regions, 
its effect is negligible. We shall therefore focus only on 
the contribution from back stress . The validity of 
this assumption is supported by the success of the recent 
single-slip theory. ^^'^^ 

Although d\j{r) is long-range, the magnitude of the 
back stress is still considerably smaller than that of 
the self-consistent internal stress rf^ when r is large com- 
pared with mean dislocation spacing. We are therefore 
interested in the contribution of d\j (r) to the stress only 
at short distances where its effect is much more pro- 
nounced. Consider a dislocation at r,, we can Taylor 



expand Kjifj) about this point, Kj{rj) ~ + {fj — 

fi) ■ Vkj + terms of higher orders. Because d\Ar) is 

symmetric while Tl^'^{f) is anti-symmetric under r i— > —r, 
the first term in the expansion vanishes. We then make 
a change of variable to the scaled coordinate r x, 
where p represents the mean total dislocations of the sys- 
tem. To second order this yields 



N 



(r-) = Y1^- j xd'ij{x)rlf{x)d'x. (47) 



Using the Fourier transform expression of dlj , the inte- 
gral in Eq. (47) can be evaluated directly using Parseval's 
theorem: 

lij^ J xd^,j{x,0)Tlf{x)d^X = J 4j{k)J^[xTlf][k]d'k 

The Fourier transform of x T^J'^ can be computed directly 
from (41): 



J^[xrlf][k] = -47rGbV^ 



{si ■ k){rhi ■ k){rhj ■ k) 



fc4 



(si • k)ipi 



(49) 



Owing to the connection d\j{x) = T fij{x), Eq. (48) be- 
comes, from (37) and (49), 



h 



A, 



d'k. 



(50) 



The vector lij is most conveniently expressed in the coordinate system of slip j. Substitution of Eq. (36) into 
Eq. (50), while projecting Sj and rhi onto {sj,rhj), gives 



(47r)^ 



-1 sin^('Afc) sm{<f>k + %) sin(30fc + 2%) 
k fc2+47rE„sin2(<^fc-^„) 



dk d(j)k 



+ m, 



2" 1 sin(0fe)sin((i6fc + %)sin(40fc-H2%) 



A 2 A: /e2 + 47rE„sin2 (</.;, -^„) 



dkd(j)k), (51) 



where Oij = {j — i)n/N is the angle between slip planes i and j. We impose a cut-off e at small k to prevent the 
logarithmic divergence due to the long-range nature of the pair correlation functions. 



r 



Under the assumption of equal interval of successive 
slip orientation, as in the previous section, we can carry 
out the above integrals very straightforwardly, giving 



lij — 



GDb 



cos{6ij)sj 



(52) 



where D = 27r^r^| lne|/A^ serves as a fitting parameter. 
The factor Xj nicely combines with p in the denomina- 
tor of Eq. (53) to make pj = Xjp. For physical rea- 
sons, we are going to replace pj with its local density 
pj{r). In the previous sections, we calculated the pair 
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correlation functions of an ensemble of spatially constant 
single-dislocation densities in thermal equilibrium. When 
the distributions of single-dislocation densities are non- 
uniform in space as is the case for systems out of equi- 
librium, the back stress response should depend on how 
much the densities vary locally. 

The final result is amazingly simple: 

rj^(r) = GDj2cosi0.,)^^^l^^^ (53) 

The above form for the back stress converges nicely to 
the single-slip theory of Groma et al.^°"^'^'^^ The cos{9ij) 
coupling between slip systems should come as no sur- 
prise. The angular dependence of the back stress must 
emerge from the symmetry of the potential. The angular 
average of tpij in Eq. (19) selects out cos{9ij) as the only 
possibility. It is interesting to note also that the same 
coupling also appears in the strain gradient theory for 
continuum crystal plasticity by Gurtin.^^"^^ 

VII. COMPARISON WITH THE EARLIER 
MULTISLIP PLASTICITY THEORY 

Recently, Yefimov et al.-'^^^^^ have proposed an exten- 
sion of their single-slip continuum plasticity theory^^'^^ 
to incorporate systems with more than one slip. In their 
theory, each slip system j contributes some amount of 
bax;k stress, given in our notation by 

r^rl = GD^^l^^j^ (54) 
to the total back stress of slip system i according to 

with slip-orientation dependent weight factor Sij acting 
as a projection matrix. For symmetry reason, three vari- 
ations were postulated: ^^'^^ 

Slj = {rhi ■ ihj){si ■ Sj) = cos^(%) (56a) 
S^j = rhi ■ [sj rhj + rhj (g) sj) ■ Si = cos(2%) (56b) 
Sfj = Si ■ Sj = cos{6ij) (56c) 

Note that the third possibility (56c) is consistent with the 
expression for the back stress we have derived in (53). 

To select among these choices, Yefimov et al. succes- 
sively used all three laws to numerically analyze the prob- 
lem of simple shearing of a crystalline strip containing 
two slip systems with impenetrable walls. The results 
of each case were compared against that from the discrete 
dislocation simulations of Shu et al."^^ The best match 
was achieved with Eq. (56b). Other choices underpre- 
dicted the amount of plastic strain. The chosen interac- 
tion law was then tested against the problem of bending 



of a single crystal strip with satisfactory agreement with 
discrete dislocation results of Cleveringa et al.^'' 

We believe that the success of their continuum theory 
in the shearing problem despite the incorrect choice of 
interaction law is due to a different reason. The amount 
of plastic strain is controlled by (i) the fitting parameter 
D and (ii) the number density of nucleation sites in the 
film. By adjusting these values, different interaction laws 
could be altered to obtain the desired fit. In their analy- 
sis, Yefimov et al. used the value of D from their previous 
single-slip theory^^ without any readjustment. There is 
no a priori reason why this value should stay unaltered. 
The density of nucleation sources in their continuum the- 
ory were chosen to match that in the discrete dislocation 
simulations. The discrepancy could also arise from dif- 
ferent ways in which the discrete dislocation theory and 
the continuum theory handle dislocation nucleation. 

In a later publication, Yefimov et al. applied their 
formalism to the problem of stress relaxation in single- 
crystal thin films on substrates subjected to thermal 
loading. Due to the difference in thermal expansion co- 
efficients between film and substrate, high tensile stresses 
can develop in the films as the temperature decreases. 
Contrary to the discrete dislocation simulations by Nicola 
et al.^^'^^ which show increasing stress built up inside a 
film with decreasing film thickness, the results from the 
continuum theory show a size-dependent hardening only 
during the early stage of cooling. Moreover, the theory 
gives identical results between some pair of slip orienta- 
tions (e.g. when the angle between the two slip planes 
6i2 is either 60° or 120°), whereas the discrete disloca- 
tion simulations and our new theory predict otherwise. 
Finally, in the previous continuum theory,^^ dislocations 
nucleate when the sum of the external stress t^'^*, the 
self-consistent long-range stress t^'^, and back stress 
exceed a certain value. From our analysis, we believe 
that, in a more correct treatment of dislocation nucle- 
ation, this back stress should be supplemented by fiow 
stress r' (Eq. (46c)) which is dominant in a nucleation 
region where plus and minus dislocations are equally pop- 
ulated. 

Applications of the current theory to the shearing 
problem and the thin film problem which shows the size- 
dependent hardening will appear shortly following this 
publication. 



VIII. DISCUSSION AND CONCLUSIONS 

We have described n^'^-order dislocation densities and 
dislocation pair correlation functions in a grand canonical 
ensemble and obtained the relationships between differ- 
ent orders of the correlation functions in the form of a 
hierarchy of integral equations. Using the Bogolyubov 
ansatz instead of the more customary Kirkwood approx- 
imation, we have closed the chain of the equations at 
second order and solved for approximate expressions of 
the pair correlation functions — valid at all distances — for 
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systems with one slip and multiple active slip systems. 
These solutions are invariant under simultaneous trans- 
formations r I— > r/ sj~p and p ^ . The transformations 
suggest that any emergent dislocation pattern should ex- 
hibit a length scale given by l/^/p as pointed out by 
Holt,^° and in agreement with the "law of similitude." 
For a complete analysis of scaling relations the reader is 
referred to Ref. 9. 

Recently Groma et al. have developed a mean- 
field variational approach to study the screening of 
dislocations/^ similar in spirit to the Debye-Hiickel the- 
ory in the study of classical plasmas. ^^"^"^ This method is 
based on approximating the system's total density ma- 
trix as a product of single-particle density matrices pi 
with the free energy given by _F = 1(H) + T'^^ Trp^ In pi. 
Although this technique provides a complimentary ap- 
proach and results in the same pair correlation expres- 
sions for a single-slip system (after some interpretation), 
its generalization to multiple-shp system is not obvious. 
In particular, one would have to supply additional cross 
couplings between different slips by hand. These cou- 
plings should automatically emerge from a complete the- 
ory. 

In Sec. VI, we have formulated transport equations for 
the total dislocation and GND densities for general mul- 
tiple slip. Interactions among dislocation pairs produce 
an additional (relatively) short-ranged "back stress" con- 
tribution to the long-range internal stress of individual 
dislocations. Most of the complexities of the correlation 
functions were integrated away, leaving only the cos{9ij) 
coupling between slip systems i and j, see Eq. (53). This 
dependence was also proposed by Gurtin in his strain gra- 
dient plasticity theory. but was abandoned by Yefi- 
mov et al.^'''^'^ We have argued in Sec. VII that this re- 
fusal was based on an unfair comparison with discrete 
dislocation simulations for the way in which dislocation 
nuncleation was treated. 

There is an important issue regarding the use of dis- 
location correlations fij for d'^ in Sec. V. The formal- 
ism developed in Sec. Ill assumes that dislocations relax 
along the directions dictated by Peach-Koehler forces. 
This implies dislocation glide, as included in the trans- 
port equations developed in Sec. VI, but also climb which 
is not considered a mechanism of plastic flow here. Math- 
ematically speaking, Eq. (15) is not the stationary state 



of Eq. (40) . Early attempts in numerically describing dis- 
location correlations in glide-only, multiple-slip systems 
failed to produce noticable patterns due to the need for 
large number of dislocations; the role of climb (or cross 
shp) was suggested to help overcome this difficulty. ^^■'^^ 
The original motivation for our approach was to find the 
orientation dependence of the back stress in the most 
straightforward way. Extracting the angular dependence 
from a climb- assisted relaxed state gave us a quick input 
to use in the glide-only multiple-slip theory. The valid- 
ity of the continuum theory will always be vindicated by 
comparisons against discrete dislocation results. 

Finally, we believe that our multiple-slip formulation 
provides a framework to address a long standing chal- 
lenge in explaining dislocation patterning. For single-slip 
systems, short-range correlations occur between two dis- 
locations except along directions normal to their glide 
plane (taken to be along y) . It has been shown that for a 
small deviation away from this "dislocation wall" direc- 
tion, an attractive parabolic potential produced by the 
correlated dislocations decays as compared with 

\y\^^ in the unscreened case.^* We have found in Sec. V, 
however, that when one or more extra slips are intro- 
duced, the effect of Debye-like screening diminishes. In 
this case, the attractive potential in fact decays like 
as if it were unscreened. This could explain the necessity 
to introduce extra slips to see the formation of walls in 
discrete dislocation simulations, -"^^"^^ unless further aided 
by climb motions. '^^'**^ The latter suggests the existence 
of a critical exponent of the attractive potential below 
which structure formation cannot occur as is the case in 
single-slip systems restricted to glide. A more detailed 
investigation of this is left for future work. 
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Note that this already implies that the analysis applies 
only to two-dimensional systems of dislocations. 
One distinction due to the choice of an emsemble type 
can be seen from the normalization condition. In the 
grand canonical ensemble, according to Eqs. (6) and 



(7), /p(")({n})4n} = EN>„^N[n:=i(,vi^ 

ni=i (N^-n V. )' where (•) denotes an average over all sta- 
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tistically equivalent ensembles. In particular, the density of 
a single dislocation of species s in a system with no external 
shear, p^^\ls) is independent of Is, thus, f p'^^^ils) d^lg = 
p^P A = (Ns). In other words p^^' = (jY,,) /A depends on 
the average number of dislocation of species s. If one were 
to carry a similar analysis in a canonical ensemble where 
the number of dislocations of each species is fixed, pi^^ 
would have to be replaced by Ns/A, where Ns is fixed. For 
higher order density, the expression becomes quite cum- 
bersome. For example, p^^J, = NgN'^/A^ for s ^ s' , while 
for s = s' , it is Ns{Ns — 1)/A'^. In thermodynamic limit 
{Ns — + oo and A ^ oo while keeping the ratio fixed) , these 
two expressions reduce to the same thing. In this sense, it 
is cleaner to work in the grand canonical ensemble. 
In the presence of an external conservative force, it can be 
shown that both (15) and (17) remain valid provided that 
an additional term representing the applied external force, 
_F(li) = — (l/fcBr)Vj^$(ri) generated by the external po- 
tential $(li) which acts on li, is added to their RHS. 
Qualitatively speaking, the original expression is nothing 
but the sum of all the Peach-Koehler interactions on the 
dislocation at li due to all other dislocations in the collec- 
tion n. 

There arc some systems where climb is typical and F is nat- 
urally small such as dislocations in vortex lattices of type 
II superconductor where the values of elastic moduli can 
be small at suitably applied magnetic field. ""^ By modifying 
the form of the interaction potential, the present analysis 
can be carried over straightforwardly. 

The form of the solution is not surprising; it suggests that 
the solution can be written as a sum of diagrams due to the 
expansion 1/(1 — a;) = 1 + x + x'^ + . . ., often encountered 
in a many-body theory. 



